/* Include files */
#include "math.h"
#include "mex.h"

/* Include to sort : http://stackoverflow.com/questions/1787996/c-library-function-to-do-sort */
#include <stdio.h>
#include <stdlib.h>


/* Checks whether inputs are matrices */
#define IS_REAL_2D_FULL_DOUBLE(P) (!mxIsComplex(P) && \
mxGetNumberOfDimensions(P) == 2 && !mxIsSparse(P) && mxIsDouble(P))
#define IS_REAL_SCALAR(P) (IS_REAL_2D_FULL_DOUBLE(P) && mxGetNumberOfElements(P) == 1)

/*************************************************/
/* SORTING */

/* Sort elements with qsort */
/*
int comp (const void * elem1, const void * elem2) 
{
	int f = *((int*)elem1);
    int s = *((int*)elem2);
    if (f > s) return  1;
    if (f < s) return -1;
    return 0;
}
*/

/* Sort index by utility */
double *base_arr;
static int compPointer (const void *a, const void *b)
{
	int aa = *((int *) a), bb = *((int *) b);
	if (base_arr[aa - 1] < base_arr[bb - 1])
		return 1;
	if (base_arr[aa - 1] == base_arr[bb - 1])
		return 0;
	if (base_arr[aa - 1] > base_arr[bb - 1])
		return -1;
	return 0;
}

/*************************************************/
/* Check Aoffer is inputing by returning in Bset */
/*
void printMatEl(double *y, double *z, int m, int n)
{	
	int A = m; int a = 0;
	int B = n; int b = 0;
	for (a = 0; a < A; a++)
		for (b = 0; b < B; b++)
		{
			z[a + A*b] = y[a + A*b];
		}
}
*/
/*************************************************/

void assignVectorWid(double *y, double *z, int M, int N, int m)
{	
	int n = 0;
	for (n = 0; n < N; n++)
	{
/*		z[n] = y[n]; */
		z[n] = y[m + M*n];
	}
}

/*************************************************/
/* MAIN */
/*************************************************/

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	/* Define output matrix */
	#define rankA	plhs[0]

	/* Define input matrices */
	#define Aoffer 	prhs[0]

	int a, b;

	/* Check the number of inputs is correct (Aoffer,Baccept,QbarA) */
	if(nrhs < 1)
		mexErrMsgTxt("Too few input arguments.");
	else if(nrhs > 1)
		mexErrMsgTxt("Too many input arguments.");

	if(!IS_REAL_2D_FULL_DOUBLE(Aoffer))
		mexErrMsgTxt("Aoffer must be real 2D full double array.");
	
	/* Dimensions and pointers for inputs */
	int A = mxGetM(Aoffer);
	int B = mxGetN(Aoffer);
	int nb = A*B;

	double *A_pt = mxGetPr(Aoffer);

	/* Output matrix */
	rankA = mxCreateDoubleMatrix(A,B,mxREAL);
	double *Rank_pt = mxGetPr(rankA);

	/* Uninitialized containers */
	mxArray *SortUtil = mxCreateDoubleMatrix(0,0,mxREAL);
	mxSetData(SortUtil, mxMalloc(sizeof(double)*1*B));
	double *Util_pt = mxGetPr(SortUtil);

	/* Generate index to use for sorting */
	int	*idx = malloc (sizeof (int) * B);
	for (b = 0; b < B; b++)
	{
		idx[b] = b + 1;
/*		printf("index %d\n", idx[b]); */
	}

	/* Assignment test */
/*	printMatEl(A_pt,Bset_pt,A,B); */

/***************************************************************/
/* Begin Sorting */
/***************************************************************/
	int ind = 0;

/* Assign pointer to vector of "offer" utilities */
	for (a = 0; a < A; a++)
	{
		assignVectorWid(A_pt,Util_pt,A,B,a);
		
	/* Sort the index based on utility values (output = idx) */
		base_arr = (double *)mxGetData(SortUtil);
		int *idxSort = idx;
		qsort (idxSort, B, sizeof(int), compPointer);

		for (b = 0; b < B; b++)
		{
			Rank_pt[a + A*b] = idxSort[b];
		}
		for (b = 0; b < B; b++)
		{
			ind = Rank_pt[a + A*b] - 1;
			if (Util_pt[ind] < 0)
			{
				Rank_pt[a + A*b] = 0;
			}
		}

/* OLD */
	/* Store the index result in the rankA matrix for each firm */
/*		for (b = 0; b < B; b++)
	    {	
	    	if (Util_pt[b] >= 0)
				Rank_pt[a + A*b] = idxSort[b];
			else
				Rank_pt[a + A*b] = 0;
	    }
*/
	}

	/* Check input */
/*  mexCallMATLAB(0,NULL,1, &Aoffer, "disp"); */

	mxDestroyArray(SortUtil);
	free(idx);

/*	return 0; */

}
